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We use N-body simulations to study the statistics of massive halos and redshift space distortions 
for theories with a standard ACDM expansion history and a galileon-type scalar field. The extra 
scalar field increases the gravitational force, leading to enhanced structure formation. We compare 
our measurements of the real space matter power spectrum and halo properties with fitting formula 
for estimating these quantities analytically. We find that a model for power spectrum, halo mass- 
function and halo bias, derived from ACDM simulations can fit the results from our simulations of 
modified gravity when erg is appropriately adjusted. We also study the redshift space distortions in 
the two point correlation function measured from these simulations, finding a difference in the ratio 
of the redshift space to real space clustering amplitude relative to standard gravity on all scales. 
We find enhanced clustering on scales r > 10 Mpc//i and increased damping of the correlation 
function for scales r < 9 Mpc//i. The boost in the clustering on large scales due to the enhanced 
gravitational forces cannot be mimicked in a standard gravity model by simply changing as- This 
result illustrates the usefulness of redshift space distortion measurements as a probe of modifications 
to General Relativity. 



The confluence of wide and deep galaxy redshift sur- 
veys with modern computing power have brought us to 
the brink of a new era for cosmology, with precision tests 
of gravity and cosmology on length scales from today's 
horizon scale to the small length scales where non-linear 
density perturbations dominate. Because galaxies and 
clusters are the tracers used to study gravity, understand- 
ing how they and their host dark matter halos form and 
evolve is crucial. The existence of structure moves galax- 
ies out of the Hubble flow, and understanding the red- 
shift space distortions, which arise from galaxy peculiar 
velocities, is an important way to extract even more in- 
formation from observations. 

Over the past two decades, great strides have been 
made in understanding the relationship between cosmo- 
logical parameters and structure formation. Large com- 
putational simulations have been performed and used to 
test and calibrate analytic approaches for understanding 
the formation of non-linear structures. Testing a theory 
of gravity that deviates from General Relativity (GR) re- 
quires checking whether the methods and results of the 
past still apply in the new model, especially beyond lin- 
ear perturbation theory. While this is true even in rela- 
tively modest alterations to gravity, such as quintessence 
dark energy models, it is particularly important when the 
new gravitational physics introduces a new "dark sec- 
tor" for gravitation that alters the gravitational force. 
Observations within the Solar System are in agreement 
with GR to great precision (for a review, see e.g. [I]). 
Hence, any new gravitational degrees of freedom must be 
suppressed on Solar System scales. There are, generally 



* markwy@oddjob.uchicago.edu 



speaking, two known ways for this to occur. 1) The effec- 
tive "charge" that responds to the new gravitational force 
is reduced by the ambient conditions in the Solar Sys- 
tem, called chameleon [2 5 J or symmetron screening |6J. 
2) The gradients that generate the new force are reduced 
through non-linear effects, which has come to be known 
as Vainshtein screening [7J [5] • Scalar fields that exhibit 
Vainshtein screening are generally called 'galileons' be- 
cause their self-interactions are determined by an inter- 
nal Galilean symmetry [5] . Because both of these screen- 
ing mechanisms are themselves inherently non-linear, it 
is important when studying them to include both den- 
sity non-linearities, as discussed above, as well as the 
non-linear structure of the modified gravity theory un- 
der consideration. This requires numerical simulations. 
Numerical simulations of galileon scalar fields will be the 
focus of this work. 

In addition to the intrinsic interest of studying what 
kinds of new gravitational strength scalar fields can exist 
in nature, there has been an important theoretical ad- 
vance in the past two years: a non-linearly complete and 
ghost-free theory of massive gravity in four dimensions 
has been found [TUHH] ■ Giving the graviton a mass adds 
a new length scale into the theory, making it possible to 
modify gravitation on long length scales in a consistent 
way. This new length scale, r c cx l/m g (where m g is 
the graviton's mass), is assumed to be of the order of 
the Hubble radius today. Solutions of this theory have 
been found that exhibit cosmological acceleration even in 
the absence of a cosmological constant [T51U5| . For the 
purposes of the present work, however, we note that this 
theory can be simplified in a decoupling limit to a the- 
ory with a ACDM background cosmology and an extra 
galileon-type scalar field [HI EH HE] that manifests Vain- 
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shtein screening [T7] . This gives additional motivation to 
our study of the cosmological effects of galileons. 

Galileon theories generally, and massive gravity in par- 
ticular, contain many of the attractive aspects of higher- 
dimensional braneworld constructions, such as the Dvali- 
Gabadadze-Porrati model [TH] and its descendants (e.g. 
[HE HO]). However, they do not require any extra dimen- 
sions and, unlike the DGP theory, are free of ghost-like 
instabilities [211 E2| . They also have the phenomcnolog- 
ical advantage that their expansion history is expected 
to be very similar to that of ACDM, whereas DGP pos- 
sessed a term linear in H in its Friedman equations that 
was difficult to reconcile with expansion history obser- 
vations, especially on the "self-accelerating" branch (see 

e.g. m)- 

The calculations we present in this work were begun 
in [231(25] in a somewhat different context. Those works 
drew their inspiration from a phenomenological version 
of a class of gravity models that arise when there are 
infinite volume extra dimensions, as occurs in the Dvali- 
Gabadadze-Porrati model and its generalizations. Those 
set-ups are generally known as "cascading gravity" mod- 
els. In [2U ES], we performed N-body simulations us- 
ing that model and characterized the non-linear power 
spectrum of the dark matter fluctuations. The models 
we considered then had a standard expansion history, 
much like the massive gravity theory; they also similarly 
increase the growth of structure on linear length scales 
while recovering GR within collapsed structures. Fortu- 
itously, the phenomenological model studied in [231 US] 
carries over nearly unchanged to the generalized galileon 
and massive gravity set-ups that are our focus in the 
present work. In brief, the results of [231 [25] were that 
semi-analytic linear perturbation theory described the 
model very well on long length scales, but that on scales 
< 10 Mpc/h, existing analytic methods for including 
non-linearities in the power spectrum failed badly. We 
also quantified, in [25], the imprint of stronger gravity 
on bulk flows, which are larger in models with an extra 
gravitational force. 

In the present work, we deepen our understanding of 
structure formation in galileon models with a standard 
ACDM expansion history by studying the power spec- 
trum, halo mass function, halo bias and the correlation 
function in redshift space measured from N-body simu- 
lations of this model. Our approach will follow the same 
pattern as work done in [26T - f2B"] for halos and [2(J1 130] for 
redshift space distortions for the f(R) model, the phe- 
nomenology of which is in some ways similar to that of 
the models we study. Our study complements and ex- 
tends similar, earlier N-body simulations performed for 
the "normal" branch of the DGP model in Ref. |31) . 

In the same spirit as [231 [25] , we do not attempt to spe- 
cialize to any particular cosmological solution for massive 
gravity or galileons (e.g. [16j[32]). Instead, we make a 
series of phenomenological assumptions designed to iso- 
late the growth-history effects of the galileon field from 
its possible modifications to the Universe's expansion his- 



tory. These assumptions are: 

• Assume exactly ACDM expansion history. 

• Compute the GR / Newtonian force as in ACDM. 

• Additionally solve the equations governing the 
galileon scalar field. 

• Assume the dynamical potential is the sum of the 
Newtonian and extra scalar field contributions. 

These assumptions mean our results will not be tied to 
any particular cosmological solution, but can be broadly 
applied to any model with galileon scalars that has an 
expansion history close to ACDM. On the other hand, 
nothing we find can be definitively associated with, for 
instance, the model of massive gravity per se, so our re- 
sults will need to be revisited as data and theoretical 
understanding of the model improve. 
Our chief results are these. 

1. The halo mass function and the linear halo bias 
at z = are altered by the increased gravitational 
force in a way degenerate with an increase in as in 
the standard ACDM model. This is in contrast 
with chameleon / f(R) theories, which generate 
a different modification to the halo mass function 
(e.g. [33,). However, since the apparent z = 
normalization of this model is itself a function of 
redshift, measurements of the halo mass function 
at different rcdshifts would break this degeneracy. 
Similar conclusions hold for the real space power 
spectrum on linear and mildly non-linear scales. 

2. Large-scale redshift space correlations are enhanced 
in these models in a way that cannot be mimicked 
within the ACDM paradigm, because velocity space 
alterations directly probe the increased long range 
gravitational force generated by the galileon field. 
In chameleon / f(R) theories, which do not have 
very long range force modifications, this effect is 
absent [25] . 

We note that the numerical approach we use is im- 
proved as compared with the N-body simulations de- 
scribed in [24 . In particular, we have replaced the phe- 
nomenological approximation used there with a more so- 
phisticated multigrid algorithm for solving the nonlinear 
equation of motion of the extra scalar field that generates 
much of the model's new physics. See Appendix [A] for 
more details. Results from this improved code first ap- 
peared in [251 133] ■ Nonetheless, the results of this much 
more computationally costly approach are in surprisingly 
good agreement with those in our previous work. 

For our numerical simulations,we will take two partic- 
ular values for the model parameter r c , which in the con- 
text of massive gravity represents the graviton's Comp- 
ton wavelength: r c = 1089 Mpc and r c = 1665 Mpc. 
These values of r c were chosen based on the results from 
[25] (see also Fig. [I]), which found (normalizing to the 
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CMB) that r c = 1089 Mpc would have a linear power 
spectrum at z — with o% — 0.92, in conflict with cur- 
rent data, while r c = 1665 Mpc would have a linear power 
spectrum at z = with erg = 0.88, which is on the bor- 
derline of being ruled out by current data. 

For our cosmology, we assume a spatially flat Universe 
with fl M = 0.24 ee fl° M , tt A = 0.76, and h = 0.73 (N.B., 
after equation 11, the notation Vl° M is used for the present 
value of SIm )■ Wherever we do not otherwise specify, we 
choose an initial amplitude for fluctuations that would 
generate erg = 0.8 in the usual GR context. We also take 
the spectral tilt n s — 0.96. Finally, wherever we do not 
write dimensionful constants explicitly, we will use units 
for which h = c = M p \ = 1. 



I. GALILEONS AND MASSIVE GRAVITY 

There are few ways to modify gravity that are not ruled 
out by Solar System tests. The simplest way to alter 
gravity phenomenologically is to add a new scalar field 
with a gravitational-strength coupling, in the spirit of the 
Brans-Dicke model. However, it is important also to look 
for modifications that have more sophisticated theoreti- 
cal motivation. Since the chief motivation for modifying 
gravity is dark energy, which only began accelerating the 
Universe's expansion rate at late times, we would like 
to find modifications to gravity that only appear at long 
length scales. 

Because of the equivalence principle, GR on its own 
cannot exhibit new behavior at long length scales. This 
is related to the fact that the graviton is massless, and 
hence has no built-in length scales other than the Planck 
scale. However, if the gravitational force were carried 
by a massive particle, the equivalence principle would 
no longer be in effect. In 1939, Fierz and Pauli demon- 
strated that a massive graviton can be defined pertur- 
batively [35]. The simple Fierz-Pauli model was found 
to be inconsistent with observations (even in the limit 
when the graviton mass is taken to zero) because of the 
so-called vDVZ (van Dam, Veltman, Zakharov) disconti- 
nuity [361 137j . This "discontinuity" arises from the phys- 
ical fact that a massive spin-2 graviton has more degrees 
of freedom than the massless one. In practice, it implies 
a disagreement between Newton's constant measured by 
gravitational lensing as compared with the gravitational 
force on massive particles. Vainshtein demonstrated that 
non-linear completions of the Fierz-Pauli evade the dis- 
continuity because the extra degrees of freedom have such 
strongly non-linear self-interactions that they decouple 
from everything else either near matter sources or as the 
graviton mass goes to zero [7], but Boulware and Deser 
showed that generic non-linear completions introduce an 
extra propagating mode that is a "ghost" (a particle with 
negative kinetic energy), implying that such theories are 
internally inconsistent [38] - However, recent work has 
shown that not every non-linear completion has a ghost. 
In particular, de Rham, Gabadadze, and Tolley (dRGT) 



found a class of completions [TUHT!?] that were shown to 
be ghost free to all orders of perturbation theory [2TJ [22] . 
As expected, the resulting theory of massive gravity prop- 
agates extra degrees of freedom beyond those present in 
GR. In particular, the graviton now has two vector com- 
ponents and one scalar component in addition to GR's 
two tensor parts. 

Although the resulting dRGT theory is both compli- 
cated and highly non-linear, at the phenomenological 
level we can gain insight by studying the theory in what 
is known as a decoupling limit. In the decoupling limit, 
we assume that the mixing among the different compo- 
nents of the graviton is small. The practical upshot of 
this is that the theory reduces to a scalar-tensor theory, 
albeit one in which the scalar field is a 'galileon' , which 
we define below [SJ [T31 [TB] . In this simpler theory, we can 
easily see the origin of the vDVZ discontinuity [3S] [37] 
and its resolution. Assuming linear theory and for wave- 
lengths small compared with the horizon scale (where the 
decoupling limit is valid), the effective Poisson equation 
for the potential, ^d yni felt by massive particles becomes 

fc 2 *dyn = -4^G(^l+0p. (1) 

The extra 1 /3 here is the manifestation of the additional 
scalar field's force. If this linear equation were exact, this 
theory would be ruled out by Solar System tests (e.g. [I]). 
However, the equation of motion for the extra scalar field, 
which we shall call ip, is strongly modified by non-linear 
effects. Making the standard and well-motivated assump- 
tion that the scalar's time derivatives are small compared 
with its spatial derivatives on sub-horizon scales, the ap- 
proximate equation of motion for tp (near flat space, i.e., 
not in an FRW-like cosmological setting) has the form 

m 

VV + |[(VV) 2 - (ViV^XVV^)] + • • • = ^T, 

(2) 

where again, in the context of massive gravity, the model 
parameter, r c = (h/c)/m g , m g is the mass of the gravi- 
ton, and r c is its associated Compton wavelength. T is 
the trace of the stress-energy tensor. The theory gener- 
ally has two more free parameters beyond r c that would 
give additional contributions in this limit. However, they 
will not be used in our study, so we simply inserted 
an ellipsis in the equation above to represent these fur- 
ther complications. Scalars with these kinds of equa- 
tions (both the simple form that appears in Eqn. [2] as 
well as its generalization) are called 'galileons' because 
their equations exhibit an analog of Galilean invariance: 
their dynamics are left the same under the replacement 
p — > p + c+b^x^; this is seen at the level of the equations 
of motion by the fact that the field is always differenti- 
ated twice. Eq. [2] is highly non-linear in the gradients 
of the scalar field, but can be solved exactly in the cases 
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of spherical symmetry and planar symmetry. A planar 
ansatz sets all of the non-linear terms to zero, whereas 
in the spherically symmetric case we find the expected 
suppression of the extra force for large densities. This 
suppression sets in at a new characteristic length scale 
r*, which is known as the Vainshtein radius; in the point 
mass case, it may be defined as 

r, = (2GA/r c 2 ) 1 / 3 = (r.r 2 ) 1 / 3 . (3) 

Roughly speaking, the extra force is unscreened at dis- 
tances larger than and is suppressed at distances 
smaller than r*. 



A. Phenomenological model 

As mentioned in the Introduction, for our numerical 
study we will make some simplifying phenomenological 
assumptions. The first of these is that we will assume 
an exactly ACDM expansion history. We do this for two 
reasons. First, we want to isolate the effect of the extra 
scalar field on the growth history for comparison with 
standard gravity. It is thus useful to keep the expan- 
sion history fixed to avoid confounding effects. Secondly, 
currently-known cosmological solutions for galileons and 
massive gravity are very close (or identical) in their ex- 
pansion history to ACDM, even in the absence of a cos- 
mological constant. Now, in addition to an expansion 
history, we must also assume some form of cosmological 
screening for the extra scalar field. That is, we expect 
the extra scalar force to be screened when the Universe's 
horizon is smaller than its own "Vainshtein" radius, and 
for the strength of the extra scalar force gradually to in- 
crease as the horizon scale grows. That is, the maximum 
strength of the force at a given time will be 

d^max — j-./ T — ^linear theory > (4) 

o±3{a) 

where B(a) is a function that depends on the cosmologi- 
cal evolution of the background. Note that this factor of 
1 /3 is the same as the one that appears in Eq. [I] when 
B = 1, we recover a scalar force whose strength is exactly 
1/3 that of GR - the maximum strength this force can 
achieve. This combination appears as the function g in 
the linearized theory of this model, and is plotted in the 
lower panel of Fig. [I] The function B(a) could be solved 
for directly in the DGP model. For the present study, we 
will follow our previous work [2H El] and adopt a form 
inspired by the DGP model: 

B(a) ^1 + 2 (Hr c f (l + JL^J . (5) 

The change in this function relative to the DGP model 
is that our function is controlled by (Hr c ) 2 , whereas the 
function in the DGP case went as Hr c . This alteration 



represents a more rapid turn-on of the scalar force rela- 
tive to DGP. Although this dependence was heuristically 
anticipated in our previous works, recent more rigorous 
attempts to study perturbations in massive gravity have 
found a function with the expected (Hr c ) 2 dependence ( 
Eq. (48) in (El) . 

Another simplifying phenomenological assumption is 
that we will drop the higher-order pieces of Eq. [2j These 
generally have the same form as the first piece, but are 
raised to higher powers. We do this mainly to make the 
numerical computations tractable, but we do not expect 
that it will make an appreciable change in the large scale 
dynamics of the theory. We can argue for this as fol- 
lows. In spherical symmetry, we can solve the equation 
exactly, with and without the extra terms (see [5]). In 
both cases, the characteristic radius that appears in the 
resulting solution has the same scaling with the mass of 
the source and r c , the graviton's Compton wavelength, 
i.e. r* cx (Mr 2 ) 1 / 3 . Hence there is no qualitatively new 
behavior or new length-scale introduced when the higher- 
order terms are included. The solutions do have different 
behavior as r — > (see e.g. |39j). but our simulations will 
not be able to resolve the length-scales on which these 
differences manifest themselves. 

We will thus solve the following system of equations: 

V 2 $at = AirGSp, (6) 

vV + ||[(vV) 2 - (v,v^)(v l v^)] = ^Sp, 

where $ is the usual Newtonian potential. Note that we 
have assumed 5p, rather than p is the source for Lp. This is 
an assumption that the phenomenological ip field we are 
studying is sourced by local overdensities, and that any 
global solutions of <p have been absorbed into the ACDM- 
like background expansion. These equations result in two 
gravitational potentials, and tp, which are combined 
into a single dynamical potential for moving particles: 

*dyn = ®N + ( 7 ) 

This is the generalization of the ^dyn that appears in 
Eq. [l] 

B. Summary of simulations 

The results reported here were compiled from a large 
number of computational runs with 512 3 particles on a 
512 3 grid. Each set of runs was comprised of simula- 
tions from z — 49 to z — (in fixed Aa = 0.0025 steps) 
performed in 4 boxes of sizes Lbox = 64, 128, 256, and 
400 Mpc/h, with three different gravitational theories: 
ordinary Newtonian gravity, plus two modified gravity 
runs assuming inverse graviton mass r c = 1089 and 1665 
Mpc. For each set, a different initialization seed was 
used to set initial conditions, then that set of initial con- 
ditions was used for each of the three different gravity 
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theories to minimize inter-run variance. We additionally 
employed a technique used in [40] that normalizes the 
overall amount of initial inhomogeneity for each initial- 
ization seed. Although these choices reduce our ability 
to directly compare our results to data, they greatly as- 
sist our comparisons with standard gravity, which is our 
focus. Some more computational details are discussed in 
Appendix [X] Our suite of simulations is as follows: 

• 8 sets of runs with outputs at only z = for all 
three models, initialized with crs(GR) = 0.8 



• 2 sets of runs with outputs at z 
and for all three models, 
cr 8 (GR) = 0.8. 



= 0.8,0.6,0.4,0.2 
initialized with 



• 1 set of GR-only runs with outputs at z = 
0.8,0.6,0.4,0.2 and 0, initialized with ct 8 (GR) = 
0.88 and cr 8 (GR) = 0.92. 

• 1 set of runs with the non-linearities in Eqs. [6] 
turned off (i.e. r c = 0) for the two galileon models. 

The initial conditions for the runs were set, as in [4"T] . 
using the Zeldovich approximation to displace particles. 
See [41] for more details. We note that our simulations 
methods are very similar to those used in [3TJ 22] ■ 



II. STRUCTURE GROWTH: LINEAR THEORY 

We can linearize the equations given in Eqs. § to give a 
simple view of how linear structure formation is altered in 
this model. In Fourier space, we can write the evolution 
of the overdensity mode 5 = Sf. as 



6" 



k 2 



a 2 H 2 



dyn ■ 



(8) 



where ' indicates the derivative with respect to e-folding 
time, din a. For linear theory, we can use 



k' 2 



*dy„ = -4ttG(1 +g)p6 k 



with g is defined by 
1 

9 = 



(9) 



(10) 



We plot g as a function of redshift in our particular cos- 
mology in Fig. [T] Using these relationships, we can find 
the combined equation 



5" 



2+*p)S' = ^n M (a)(l + g)6 (11) 



where, for easy reference, 
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FIG. 1. Results from linear theory. In the top panel, we 
show the z — amplitude of the dark matter linear power 
spectrum that would be inferred by standard methods for 
different values of r c . This measure, Apparent <7g, is defined 
in the text and in the plot title where Dgr is the growth 
factor in standard gravity. In the lower panel, we plot the 
ratio of the strength (in linear theory) of the extra scalar 
force as compared with the Newtonian force; this appears as 
g = l/(3B(a)) in the text. Although the extra force has in 
principle a maximum strength of 1/3 the Newtonian force, the 
phenomenological model we adopt suppresses this under the 
assumption that the background density of space will modu- 
late the strength of the extra scalar's force, as was found in 
the DGP model. 



For a cosmology with only matter and A, we have 

H' 3 

~H ~ ~ ~ 



n° M a- 3 



+ (i - n° M ) 



(12) 



For our linear theory solutions, we take initial condi- 
tions from CAMB @3] at z = 50. At this redshift, the 
extra scalar has no effect on cosmology, so we are justi- 
fied in using a standard gravity-based code for generating 
our initial conditions. We then use Eq. [TT]to numerically 
evolve the Fourier density modes from z = 50 to z = 
for different values of r c . We plot the results in Fig. [l] 

The results plotted in Fig. [T] are presented in a some- 
what unusual format, which we will explain. As we have 
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stated before, the extra scalar force only starts to in- 
fluence growth history at relatively late times, and then 
grows in its influence as time goes on. Hence, it is useful 
to quantify the deviation of the power spectrum gener- 
ated including the extra scalar from the power spectrum 
that would have been generated by GR. We do this by 
means of a measure we call Apparent z — a%. Re- 
call that a 8 is a commonly used quantity that encodes 
the normalization of the power spectrum by measuring 
the matter fluctuation within 8 Mpc/h spheres. In gen- 
eral, we can define the variance of the linear density field 
within spheres of radius R, <jr via 



2 



1 

2^2 



k 2 P{k, z = 0)W|(fc)dfc . (13) 



To get of, we take W s {k) = 3ji(kR 8 )/kR 8 , with R s = 8 
Mpc/h and j\ a spherical Bessel function. 

Observational cosmology involves making measure- 
ments of the power spectrum of density perturbations 
at many redshifts. For comparison purposes, these mea- 
surements are commonly extrapolated to z — using 
linear theory and reported as measurements of og. It is 
in this sense that the amplitude of the CMB power spec- 
trum gives a measurement of Og. In Fig. [T] we make use 
of a similar algorithm: 

• For a set of Zj between 5 and 0, we evolve the power 
spectrum from z; n j t = 50 to z — Zj using the mod- 



ified gravity equations, Eq. 11 



• Beginning at z = Zj , we take the modified gravity- 
generated power spectrum and then evolve it to 
z = using the growth factor, -DgRj i n standard 
gravity, essentially setting g = in Eq. [Tl] 

• We evaluate o s using the power spectrum generated 
by this procedure, and record it as the Apparent 
a s (zj). 

These linear theory results give a simple qualitative 
picture of how this model modifies the growth of struc- 
ture: 

1. Growth is the same as in GR when the Universe's 
horizon scale is smaller than its Vainshtein radius. 

2. As the Universe expands, its horizon grows larger 
than its Vainshtein radius, allowing the extra scalar 
to begin to accelerate the growth of structure (see 
the second panel of Fig. [TJ 

3. As the extra scalar force operates, structure grows 
faster than it would in GR, leading to more large 
scale structure and a larger Apparent <jg. 



III. RESULTS: HALO PROPERTIES 

In this section we study the halo abundance and linear 
bias, as inferred from simulations, and how standard fit- 
ting prescriptions calibrated in ACDM simulations can 



match the galileon modified gravity simulated results. 
The techniques we employ to extract the halo catalog, 
mass-function and bias are the same as those presented 
in [55] in the context of f(R) models. 



A. Halo Abundance 

We detect halos using a spherical overdensity (SO) 
halo finder, and define halo masses within an overden- 
sity A = 200 with respect to the mean background den- 
sity p m . The halo mass function is obtained dividing the 
number of halos in a given mass bin by the comoving 
volume of the simulation box and the mass bin size. To 
be conservative, we keep only halos with more than 6400 
particles in each box and combine the mass functions 
from all boxes and runs at z — 0. 

We compare our measured mass function with the fit- 
ting formula of Tinker et al. [44], obtained from high- 
resolution ACDM simulations, and given by 



n M = 



dn(M, z) 
dlogAf 



diner 



M dlogM 



(14) 



where a(M) = an is the variance of the linear density 
field for a mass M — 47ri? 3 / o TO /3 contained in a sphere of 
radius R at the mean background density (see Eqn 
and 



131, 



/(*) = A 





— a 








hi 







-c/rr 2 



(15) 



For A = 200 we set parameter values A = 0.186, a 
1 .47. 6 = 2.57 and c= 1.19. 

In Fig. [2j we show the mass function as a function 
of halo mass for simulations and predictions. In the 
panels on the left, the simulation-derived mass function 
is plotted as triangles for the ACDM simulations and 
squares for galileon modified gravity simulations with 
r c = 1665 Mpc (top) and r c = 1089 Mpc (bottom). Mean 
values and error bars are derived from volume-weighted 
bootstrap samples in order to reduce sample variance, 
similarly to [37J 12H]- In the right hand panels, we show 
the relative change AnAf/nju = {n^/n 1 ^ 13 ^ 1 — 1) mea- 
sured from the r c = 1665 Mpc (triangles) and r c = 
1089 Mpc (squares) simulations. Here we have neglected 
the last measured point shown on the left panels, corre- 
sponding to the most massive halos, because the number 
of halos in this bin is of order unity, making any measure 
of relative differences meaningless. We also show the pre- 
dicted percent difference between the Tinker et al. [H] 
fit for ACDM with a$ = 0.88 and 0.92 relative to that 
for a$ = 0.8. 

We find that the mass function in the modified grav- 
ity scenario can be fit with a ACDM mass function with 
higher erg. This degeneracy prevents low redshift clus- 
ter abundance measurements from being able to distin- 
guish galileon modified gravity from ACDM. In principle 
this degeneracy can be broken by cluster observations at 
higher redshifts. As we go back in time, the extra force 
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FIG. 2. (Left): Halo mass function at z — measured from ACDM (red triangles) and galileon modified gravity simulations 
(blue squares) with r c — 1665 Mpc (upper panel) and 1089 Mpc (lower panel). The ACDM predictions of Tinker et al. [44] for 
<7 S = 0.8 (solid curves) and as = 0.88,0.92 (dashed curves) are also plotted. (Right): Percent difference for the mass function 
measured from modified gravity simulations relative to a ACDM simulation with as — 0.8 (triangles and squares), and Tinker 
predictions relative to as = 0.8 (dashed lines). 



is weaker and the cluster measurements tend to agree 
with those of a true ACDM universe (in our case with 
as = 0.8). From Fig. [I] we see for instance that the ap- 
parent erg for the r c = 1089 Mpc case changes from 0.92 
at z = to 0.86 at z = 2. For higher values of r c , this 
change is smaller and it becomes harder to distinguish 
galileons from ACDM. 



B. Halo Bias 

Halo bias characterizes the clustering of galaxies with 
respect to the underlying mass distribution, and is as- 
sumed to be scale independent on large scales which are 
still in the linear regime. For halos of a given mass M, 
we compute the halo bias from the simulations using the 
definition b(k,M) = Phm{k) / P mm {k) , where P mm (k) is 
the dark matter power spectrum and Phm (k) is the halo- 
mass cross spectrum. This choice allows us to partially 
reduce the shot noise that would result if we used the 
halo-halo auto-spectrum Phh{k). To obtain the linear 
halo bias, we fit a polynomial to b(k, M) using the first 
10 values of k and extrapolate to the lowest k value, i.e. 
b L (M) = b(k = k min ,M). 

Similarly to the mass-function results, we then com- 
pare our measured bias with the fitting formula of Tinker 
et al. [IS] , calibrated from high-resolution ACDM simula- 
tions as 



b L (M) = 1 — A- 



- Si 



Bv b + Gv c 



(16) 



where v(M) = S c /a(M) and S c — 1.686. We fix values 
for parameters A, a, B, b, C and c appropriate for our 
SO mass definition of A = 200 05J. 

In Fig. [3j we show the linear halo bias for the same 
cases shown in Fig. [2] The bias measurements are nois- 
ier than those of the mass-function, mainly due to the 
shot noise in the halo-mass cross-spectrum. The per- 
cent difference Ab L /b L = (b r £ /&£ CDM - 1) is less pro- 
nounced than the corresponding differences in the abun- 
dance. Nonetheless, it is consistent with a ACDM model 
with a higher value of <7g. 

Altogether, the halo properties at low z show that it 
should be hard to break the degeneracy between galileon 
modified gravity and ACDM models with larger values 
of cr 8 unless one can make measurements at various red- 
shifts. Even though it is beyond the scope of this work, 
it should be possible to develop improved methods that 
break universality in order to include the effects of mod- 
ified gravity more accurately in the halo properties. In 
the remaining of the paper we investigate whether a more 
unique signature of galileon models can be seen in the 
two-point statistics of galaxies and in their redshift-space 
distortions. 



IV. RESULTS: CLUSTERING IN REAL AND 
REDSHIFT SPACE 

In this section we present the power spectrum mea- 
sured from the simulations in real space (Section IV A) 



and review the linear perturbation theory of redshift 




FIG. 3. Linear halo bias at z = for the same cases shown in Fig. [2] The results for the halo bias are noisier than those for 
the mass function, but similarly display the degeneracy between galileon modified gravity and ACDM with larger ag at single 
redshifts. 



space distortions (Section IV B |. The method used to es- HALOFIT prediction at erg = 0.8. 



timate the correlation function in real and redshift space 
is outlined in Appendix [SJ Our results showing the red- 
shift space clustering signal in galileon models compared 
to ACDM are presented in Section [IV C| 



A. Power spectra in real space 

In order to estimate the power spectra, for each box of 
a given size and cosmological model, we define a density 
field in the grid of 512 3 points using the Cloud-In-Cell 
(CIC) method. The density field is Fourier transformed 
and the power spectrum of each box is estimated av- 
eraging the band-power in each Fourier mode bin. We 
then combine the various boxes using the same volume- 
weighted bootstrap averaging procedure used for the halo 
mass-function and bias. 

In Fig. [4] we show the power spectra for the same cases 
displayed in Figs. [2] and [3j On the left panels we can 
verify that the ACDM simulations shown with red trian- 
gles are well described by the HALOFIT fitting formula 
(46] . shown as a solid black line; for completeness the 
dotted line shows the linear power spectrum in this case. 
The results for the galileon simulations are shown as blue 
squares for r c = 1665 Mpc (top) and r c = 1089 Mpc 
(bottom). We also display in dashed lines the results 
from HALOFIT with different values of a a . 

On the right panels we show percent differences 
AP/P = (p'v/pACDM _ for thege cages ^ The trian _ 

gles (squares) compare the r c = 1665 (1089) Mpc simu- 
lations to the ACDM simulation. The dashed lines com- 
pare the HALOFIT prediction at er 8 = 0.88 (0.92) to the 



Linear theory predicts a scale-independent deviation 
of - 20% and - 32% for r c = 1665 Mpc and 1089 Mpc 
respectively, relative to ACDM. This is in fact observed 
in the simulations on the largest linear scales. For mildly 
nonlinear scales, the HALOFIT prescription with higher 
<78 provides a close description of the departures seen 
in the galileon simulations relative to those of ACDM 
simulations. The departures are not well captured for 
scales that are more deeply inside the non-linear regime, 
due to the Vainshtein mechanism. On these scales, how- 
ever, baryonic physical processes (not considered in our 
results) arc in effect and represent an extra source of de- 
generacy. 



Therefore, our results from the power spectra in real 
space at z = indicate that on linear and mildly non- 
linear scales there is a degeneracy between galileons mod- 
ified gravity and ACDM with higher erg, similarly to that 
seen for the halo mass-function and linear bias. Again, 
this degeneracy is redshift dependent and measurements 
of the power spectrum at distinct redshifts could in prin- 
ciple isolate features of galileon models. On deeply non- 
linear scales, the Vainshtein mechanism brings about 
unique features of galileon models at a single redshift, but 
processes due to the baryonic physics become more rele- 
vant, likely representing an even more important source 
of degeneracy, unless such complex processes are well 
characterized. 
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FIG. 4. Power spectrum at z — for the same cases shown in Figs. [2] and [3] The left panels show the simulated power spectra 
for the ACDM simulation (triangles) and for the galileon simulations (squares) with r c = 1665 Mpc (top) and r c = 1089 Mpc 
(bottom). Also shown are the linear power spectrum for ACDM with erg = 0.8 and the HALOFIT power spectra for erg = 0.8 
(solid line) 0.88 (dashed top) and 0.92 (dashed bottom). The right panels show relative deviations from these cases. On deeply 
non-linear scales the Vainshtein screening is effective and it is not possible to fit the galileon simulations using standard ACDM 
fits with higher values of erg. 



B. Modeling redshift space distortions 

The growth rate of structure in the Universe may 
be determined through the observed anisotropy of the 
galaxy clustering in redshift space, caused by the line of 
sight component of the galaxies peculiar velocities. Red- 
shift space effects alter the appearance of the clustering 
of matter on all scales, and together with nonlinear evolu- 
tion and bias, give rise to the measured anisotropic power 
spectrum or correlation function which is different from 
the simple predictions of linear perturbation theory. The 
comoving distance to a galaxy, s, differs from its true 
distance, x, due to its peculiar velocity, v(x) (i.e. an ad- 
ditional velocity to the Hubble flow) . The mapping from 
redshift space to real space is given by 



u z z, 



(17) 



where u z = v-z/(aiJ) and H(a) is the Hubble parameter. 
This assumes that the distortions take place along the 
line of sight denoted by z. Note this is the plane parallel 
approximation which we adopt in this paper. 

On large scales, coherent infall into overdense regions 
distorts clustering statistics, causing the correlation func- 
tion to appear squashed along the line of sight [see|37J for 
a review of redshift space distortions] . For growing per- 
turbations in the linear perturbation regime, the overall 
effect of redshift space distortions is to enhance the clus- 
tering amplitude. This can be seen as an enhancement of 
the power spectrum in redshift space, P s (k), compared 



to that in real space, P r {k). This effect was first ana- 
lyzed by [15] in linear perturbation theory and can be 
approximated by 



(18) 



where [i is the cosine of the angle between the wavevector, 
k, and the line of sight. The variable j3 is 



n 1 dlnZ? / 

f3 = = - 

b dlna b 



(19) 



where / is referred to as the linear growth rate and b is the 
linear bias. In this section we restrict the analysis to large 
scales and assume a constant bias b = \/ £,hh{ r ) / £,m{r) 
where £hh(m) is the two point correlation function for the 



halos (dark matter). The 'Kaiser formula' (Eq. 18) re- 



lates the overdensity in redshift space to the correspond- 
ing value in real space and is the result of several approxi- 
mations, e.g. that the velocity and density perturbations 
satisfy the linear continuity equation. All of these as- 
sumptions are valid on scales that are well described by 
linear perturbation theory and will break down on differ- 
ent scales as the density fluctuations grow [see e.g. l49Tf5l| 
for more details]. 

Rather than use the full 2D power spectrum, P(k,fj,), 
it is common to decompose the matter power spectrum 
in redshift space into multipole moments using Legendre 
polynomials, Li(n), [see e.g. 27] 



P(k,fi)=Y,Pi(k)Li(fi), 



(20) 
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where the summation is over the order, I, of the mul- 
tipole. The anisotropy in P(k) is symmetric in fi, as 
P(k, n) = P(k, — /Lt), so only even values of I are summed 
over. Each multipolc moment is given by 



Pi s (k) = 



21 + 1 



P(k,ii)Li(ji)dfj,, 



(21) 



where the hrst two non-zero moments have Legendre 
polynomials, Lq(ji) — 1 and L 2 (/i) = (3fi 2 — l)/2. Using 



the linear model in Eq. 18 the first multipole moment is 
given by 



1 



Po(k) = P m (k)(l + -/3+-p 2 ) 



(22) 



where P m (k) denotes the real space matter power spec- 
trum. Note we have omitted the superscript s here for 
clarity. 

The corresponding equation in configuration space can 



be obtained by Fourier transforming Eq. 22 giving the 
corresponding relation between the monopole of the cor- 
relation function in redshift space to real space [57] 



(23) 



where £(r) is the real space correlation function. The 
above equations describe the boost in the clustering sig- 
nal in redshift space on large scales where linear pertur- 
bation theory is valid. To go beyond linear theory and 
deal with small scale velocities requires a model for the 
velocity field and all the density velocity correlations. On 
small scales, randomized velocities associated with the 
motion of galaxies inside virialized structures reduce the 
power. The dense central regions of galaxy clusters ap- 
pear elongated along the line of sight in redshift space, 
which produces the 'fingers of God' (FOG) effect seen in 
redshift survey plots [52]. This FOG effect can be de- 
scribed by convolving the correlation function £(sj_, sm), 
where sh is the distance separation along the line of sight 
and s± is the perpendicular separation, with the distri- 
bution function of random pairwise velocities, f(u) [53], 



£( S -L, S ||) = 



du/(u)£(sx,S|| - u) 



(24) 



where f(u) can have an exponential or a Gaussian form 
such as 



/(«) = 



1 



V 7 ^ 



?exp 



u 

2^2 



(25) 



where <r v is the pairwise peculiar velocity dispersion. 
This model has been used to fit to results from both simu- 
lations and observations [see, for example 14TH [5TJ1 f54l!57| . 
Recently there have been many models which improve 
on this description of redshift space distortions in the 
nonlinear regime [49 1 [58H61] . In this work our main con- 
cern is to quantify the relative difference in the correla- 
tion function measured from the simulations in redshift 



and real space in ACDM compared to galileon models. 
We restrict our study to large scales and shall compare 
the measured ratio between the redshift and real space 
correlation function with the linear perturbation theory 
predictions in each cosmology. 

Linear theory results for various values of r c 
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FIG. 5. The linear growth rate - both weighted by <Jg,(z) 
and alone - as a function of redshift, z, for ACDM, the 
r c — 1089 Mpc and r c — 1665 Mpc galileon models, plus 
a comparison model with r c = 4000 Mpc, are shown as a 
black dashed line, a solid black line, a dashed blue line, and 
a dash-dotted red line, respectively. The grey band around 
the ACDM line represents the expected precision of future 
surveys, which hope to achieve ±2% accuracy in their mea- 
surement of the growth rate. 



C. Redshift space distortions from galileons 

Measuring the anisotropic distortions in the galaxy 
clustering pattern in redshift space constrains the pa- 
rameter /3 = f /b, where b is the galaxy bias factor and / 
is the logarithmic derivative of the linear growth rate of 
structure, which is scale independent in the case of gen- 
eral relativity. In the lower panel of Fig. [5] we show the 
linear growth rate as a function of redshift for ACDM 
and the r c = 1089 Mpc and the r c — 1665 Mpc galileon 
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models as a dashed black line, a solid black line and a 
dashed blue line respectively. A model with r c = 4000 
Mpc is also plotted as dot-dashed red line for compari- 
son. In the upper panel we plot the linear growth rate 
weighted by cs(z) for each model. The grey shaded re- 
gion around the ACDM result represents the expected 
precision of a DETF [52] Stage IV galaxy redshift survey 
such as the ESA's EUCLID mission [55], WFIRST [55] 
or the ground-based dark energy experiment, BigBOSS 
[65], which aim to achieve ~ 2% accuracy in their mea- 
surement of the growth rate. The relative difference in 
the linear growth rate between the r c — 1089 Mpc model 
and ACDM varies from 16% at z = to 12% at z = 1, 
while for r c — 1665 Mpc the difference is 12% and 8% at 
redshifts z = and z = 1 respectively. It is clear that 
if future galaxy redshift surveys can measure the growth 
rate to within 2%, they will be able to place significant 
constraints on currently allowed galileon modified gravity 
models. 

In Fig. [6] we show the two point correlation function 
in real and redshift space measured from the simulations 
at z = for ACDM and the r c = 1089 Mpc (r c = 1665 
Mpc) model in the left (right) panels. In the upper panels 
open symbols represent the correlation function in real 
space while closed symbols are used for the redshift space 
function £(s). It is clear from this figure that both the 
real and redshift space correlation function amplitude are 
increased in the modified gravity models compared to 
ACDM. We plot the relative difference in the ratio of the 
correlation function measured in redshift space to that in 
real space, in the r c — 1089 Mpc (r c = 1665 Mpc) model 
compared to ACDM in the lower left (right) panel. Here 

A(£s/£r) = (&,/Cr)r c /(6/£r)ACDM - 1 where (£ a /£r)r IS 

the ratio of the monopole of the redshift to real space 
correlation function for the modified gravity model. 

From Fig. [6] we find an increase in the clustering signal 
in redshift space on large scales in the modified gravity 
model compared to ACDM and an increase in the small 
scale damping due to incoherent random velocities. This 
increased clustering signal at r > 10 Mpc/h in the mod- 
ified gravity models is due to increased bulk flows on 
large scales. On small scales the enhanced forces in the 
modified gravity model create a larger velocity dispersion 
which gives rise to increase damping compared to ACDM 
on scales r < 10 Mpc/h. These results agree with a sim- 
ilar study of redshift space distortions in f(R) modified 
gravity carried out by (29) . The Kaiser linear theory pre- 
diction for this ratio using the appropriate growth rate 
for each model is shown in the lower panels in Fig. [6] as a 
dotted grey line. The measured difference in the redshift 
to real space ratio between ACDM and galileon mod- 
els agree with linear theory predictions on scales r > 15 
Mpc/h for the r c = 1089 Mpc model and r > 8 Mpc/h 
for r c — 1665 Mpc. 

Both the upper and lower panels in Fig. [6] combine 
measurements from the L box = 400 Mpc/h and L box = 
256 Mpc/h simulations. These measurements represent 
the average over 8 realizations and the errors plotted rep- 



resent the scatter amongst these 8 different simulations 
for each cosmological model. We plot the ratio as in 
the lower panels in Fig. [6] to remove any sample vari- 
ance in the measurements which arise from sampling a 
finite number of large scale modes in a finite simulation 
volume. 

In Fig. [7] we plot the relative difference, as in the 
lower panel in Fig. [6j for the r c = 1089 Mpc (left) and 
r c = 1665 Mpc (right) model at z = 0.4, z = 0.6 and 
z = 0.8 in the lower, middle and upper panels respec- 
tively. The relative difference in the ratios increases with 
redshift from e.g. 4% at z — to 6% at z = 0.8 for the 
r c = 1089 Mpc model and agrees with the linear per- 
turbation theory predictions (grey dotted line) on larger 
scales compared to z = 0. The enhanced forces due to 
the galileon scalar field decrease with increasing redshift 
and as a result the relative increase in the small scale 
damping decreases. The differences between the modi- 
fied gravity models considered here and standard gravity 
decrease with increasing redshift and so it may be coun- 
terintuitive that the ratios shown in Fig. [7] increase with 
redshift on large scales. The reason for this can be found 
by examining the linear perturbation theory function for 
the ratio £ s /£ r = 1 + 2/3/ + 1/5/ 2 for the dark matter. 
For a given cosmological model, this ratio increases with 
increasing / (increasing redshift). Even though the rela- 
tive difference in / between a modified gravity cosmology 
and ACDM decreases from z = to z = 0.8, the relative 
increase in the ratio will be larger at z = 0.8 compared 
to z = 0, as can been seen from Fig. [7] The errors plot- 
ted here represent the Jackknife errors on the mean as 
outlined in Appendix [B] 

From Figs. [2] and [3] it is clear that both the halo mass 
function and linear bias, calculated using the Tinker fit- 
ting formulae [44 , 45] can match the measurements from 
the r c — 1089 Mpc and the r c — 1665 Mpc simulations if 
we assume a ACDM cosmology with a higher value for erg. 
Moreover, from Fig. [4] the non-linear real-space power 
spectrum calculated using the HALOFIT fitting formula 
|46j can also match the modified gravity at mildly non- 
linear scales. We test if this degeneracy also occurs in 
the measurements of the correlation function in redshift 
space by running two additional ACDM simulations us- 
ing the same O m but with erg = 0.92 and ag = 0.88 at 
2 = 0. In Fig. [8] we plot the same ratio of £ s /£ r for 
the r c = 1089 Mpc (r c = 1665 Mpc) model compared 
to ACDM in the top (bottom) panel as a red solid line 
with error bars, as shown in Fig. [6j The relative differ- 
ence in £ s /£ r measured from a ACDM simulation with 
(j g = 0.92 (<rg = 0.88) compared to a ACDM simulation 
with erg = 0.8 is shown as a black dashed line in the top 
(bottom) panel. The grey shaded regions represent the 
jackknife errors on the mean. 

The difference in the monopole to real space corre- 
lation function between ACDM cosmologies which have 
different power spectrum amplitudes at z = has a dif- 
ferent signature on large scales compared to the mea- 
surements from our modified gravity simulations. The 
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FIG. 6. Upper panels: The dark matter two point correlation function measured in real and redshift space at z — are shown as 
red circles and blue squares for ACDM and the r c = 1089 Mpc (r c = 1665 Mpc) model in the left (right) panel. Open symbols 
denote the correlation function measured in real space while closed symbols represent the correlation function measured in 
redshift space. Lower panels: The relative difference in the ratio of the correlation function measured in redshift space to that 
in real space, in the r c — 1089 Mpc (r c = 1665 Mpc) model compared to ACDM is shown in the left (right) panel. The dotted 
grey line in both the right and left panel shows the Kaiser linear theory prediction for this ratio using the appropriate growth 
rate for each model. 




r [Mpc/h] r [Mpc/h] 



FIG. 7. The relative difference in the ratio of the correlation function measured in redshift space to that in real space, in the 
r c = 1089 Mpc (r c = 1665 Mpc) model compared to ACDM is shown in the left (right) panel. The upper, middle and lower 
panels show this ratio measured from the simulatons at z = 0.8, z = 0.6 and z = 0.4 respectively with Jackknife errors on the 
mean. The dotted grey line in both the right and left panels shows the Kaiser linear theory prediction for this ratio using the 
appropriate growth rate for each model at the redshift indicated by the legend. 
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The dotted black line represents the predictions of 
Eq. 23 where (3 — f/b and we use the linear theory 
value for the growth rate in each cosmology at each red- 
shift. We take the line ar bias to be th e best fit value 
for the ratio b(M) = ^hh (r, M) / £ m (r) over the range 
r = 20 — 50Mpc//i, where £,hh{.f, M) is the halo correla- 
tion function in real space and £ m (r) is the dark matter 
correlation function. We find that the bias for ACDM 
halos in the mass range > 10 13 M & /h varies from b m 2 
at z — to b w 3.9 at z = 0.8 compared to the r c = 1089 
Mpc model which has b = 1.5 (b = 2.9) at z = 0(z = 0.8). 
The measured relative ratios in Fig. [9] are consistent with 
the linear perturbation theory prediction at all redshifts 
within the error bars with a 4% (5%) difference between 
the r c = 1089 Mpc model and ACDM at z = (z = 0.8). 
A similar analysis for coupled dark energy cosmologies 
was recently carried out by |68j where they also exam- 
ined the two point correlation function of halos in redshift 
and real space at several redshifts. 



FIG. 8. The relative difference in the ratio of the correlation 
function measured in redshift space to that in real space, in 
the r c = 1089 Mpc (r c = 1665 Mpc) model compared to 
ACDM is shown in the top (bottom) panel as a red solid 
line with error bars as in Fig. [6] The difference in the ratio 
£ s /£r measured from a ACDM simulation with as = 0.92 
(as = 0.88) to a ACDM simulation with as = 0.8 is shown as 
a black dashed line in the left (right) panel. The grey shaded 
regions represent the jackknife errors on the mean. 



ratio on large scales is consistent with unity and is dis- 
tinguishable from the relative increase in this ratio mea- 
sured in the modified gravity simulations compared to 
ACDM, while on small scales we measure a similiar in- 
crease in the damping signal in the ACDM simulation 
with as = 0.92 (erg = 0.88) compared to the simulation 
with as = 0.8. In linear perturbation theory the velocity, 
v oc fa 8 and in this case the growth rate, /, is the same 
for all the ACDM simulations. This implies that the ve- 
locities for the us = 0.88 (0.92) simulation will be higher 
then in the as = 0.8 case and explains the decrease in 
the ratio on small scales shown as a black dashed line 
m Fig. [8j These results agree with similar studies of the 
correlation function in real and redshift space for ACDM 
cosmologies varying the parameter combination (f2 m , as) 
which were carried by 66. 167). 

In Fig. [9] we plot the ratio of £ s /£ r measured using 
all halos with masses > 10 13 M Q //i in the r c = 1089 
Mpc (r c = 1665 Mpc) simulation compared to ACDM 
in the left (right) panels as a green solid line. The mean 
value plotted in this figure is from two simulations and 
the gray shaded region represents the Jackknife errors on 
this mean. The resolution of our simulations limits us 
to using only the 400 Mpc/h simulation boxes for this 
measurement and so we cannot probe the difference in 
the nonlinear damping signals for these halo mass ranges 
at r < 10 Mpc/h. 



V. CONCLUSIONS 

We have studied linear and non-linear structure for- 
mation in a class of modified gravity models with a 
ACDM-like expansion history and a galileon scalar field 
that is screened in regions of high density via the Vain- 
shtein mechanism. Our primary results have been de- 
rived from a large suite of N-body simulations of this 
model, where the effects of the extra scalar field have 
been captured through direct numerical solution of its 
non-linear equations of motion. We have found, as ex- 
pected, that large scale structure is enhanced in this 
model relative to ordinary general relativity. This en- 
hancement appears through a significant increase in the 
number of large halos formed in this model relative to 
standard gravity, with the enhancement matching that 
anticipated from linear theory: the modified gravity mod- 
els lead to a late-time power spectrum normalization of 
ag = 0.92 (0.88) for the values of r c that we studied, 
r c = 1089 (1665) Mpc. Hence, for a single redshift, the 
effect of the galileon field is degenerate with an increase in 
the primordial amplitude of fluctuations. However, as we 
emphasized in Fig. [l] the apparent z — normalization 
of the power spectrum is itself a function of redshift in 
this theory, so measurements of the halo mass function at 
different redshifts would break this degeneracy. This de- 
generacy also appears in the real space power spectrum 
even at mildly non-linear scales. On deeply non-linear 
scales the Vainshtein screening provides a unique feature 
of galileon models, but is likely degenerate with the much 
more complex baryonic physics that becomes more rele- 
vant on these scales. 

Next, we examined the impact of redshift space distor- 
tions on the measured two point correlation function in 
our simulations. We find deviations in the ratio of the 
redshift to real space correlation function £ s /£ r in the 
modified gravity models that are clearly distinguishable 
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FIG. 9. The ratio of £ s /£r measured using all halos with masses > 10 13 Mq /h in the r c — 1089 Mpc (r c = 1665 Mpc) model 
compared to ACDM are shown in the left (right) panels as a green solid line. The results at z = 0, z — 0.4 and z = 0.6 are 
shown in the bottom, middle and top panels respectively. The grey shaded region show the Jackknife errors on the mean. 



from standard gravity. On large scales the redshift space 
correlation function is sensitive to the growth rate, which 
for the modified gravity models we consider here, can be 
dramatically different from the growth rate in standard 
gravity. For the values of r c that we studied, we found 
deviations in the dark matter clustering signal on large 
length scales at the level of 10% and a difference of 4-5% 
at z = — 0.8 for halos > 1O 13 M //i which is potentially 
large enough to be seen or ruled out by future galaxy 
redshift surveys. As we show in Fig.|HJ the redshift space 
distortion signal on large scales cannot be mimicked by 
increasing the initial amplitude of fluctuations and using 
standard gravity. Going beyond what we can model in 
linear theory, we also find that the enhanced gravitational 
force gives rise to a diminution of the clustering signal on 
small (< 9 Mpc/h) length scales beyond that seen in GR, 
as the enhanced non-linearities more efficiently wash out 
the correlated motions found on large length scales. 

Let us compare our results to those found in f(R) grav- 
ity. In contrast with our findings, the halo mass functions 
for f(R) models are not well approximated by simple 
changes in erg, even for a single redshift. Furthermore, 
since the extra scalar field that appears in f(R) is always 
massive, it necessarily has a range limited by Yukawa 
suppression. Hence, f(R) models cannot generate the 
long-range deviations from GR that are present in the 
redshift space correlations we have studied in galileon 
models. This is an important illustration of how f(R) 
gravity and galileon models are quite distinct in their 
phenomenology. 



Despite the presence of Vainshtein screening in regions 
of high density, our findings suggest that linear perturba- 
tion theory is a good guide to understanding how large 
scale structure is modified in models with a galileon scalar 
field, at least for the values of r c we have studied. If, as 
data improve, the limits on r c are pushed upwards to- 
wards today's Hubble scale, the methods we employ will 
have to improve accordingly to make further progress. 
On the theoretical side, it will be necessary to under- 
stand how to improve the phenomenological model de- 
scribed in §1 A| to represent better the interplay between 
the cosmological and local excitations of the scalar mode 
of the graviton. On the computational side, our simula- 
tion methods will have to improve in order to capture a 
greater range of scales so that we can find and quantify 
the scalar's effects as they become more subtle, as they 
will with larger r c . Since r c ~ c/Hq is the value we ex- 
pect if today's cosmological acceleration is generated by 
a massive graviton, we will only be able to use the growth 
of structure to constrain the theory definitively once all 
of these improvements are made. 

Note added: While this manuscript was in the final 
stages of preparation, Ref. jSH] appeared, describing a 
new code that can solve galileon-type equations on an 
adaptively refined mesh. However, they study the model 
of self-accelerating DGP, which has a different expansion 
history from ACDM and a repulsive, rather than an at- 
tractive, extra scalar force. Hence their results are phys- 
ically distinct from the phenomenological model studied 
in this paper. 
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Appendix A: Numerical details 

For our N-body simulations, we improved the code first 
reported on in 24., solving a very similar set of equa- 
tions as those studied in [42]. The code is written in 
FORTRAN and is largely based on the publicly available 
PM-Code [JT] , a particle-mesh N-body code that employs 
fast fourier transforms to solve the Newtonian Poisson 
equation and cloud-in-cell grid assignment to interpolate 
discrete particle positions onto the density grid. The pri- 
mary addition we have made to the public code beyond 
those reported in |24j is the inclusion of a multigrid relax- 
ation subroutine for solving the non-linear equation for 
the extra scalar field, Eqs. [6] This subroutine is heavily 
adapted from [7D]. We have also included threaded par- 
allelization through OpenMP. The numerical method we 
utilize is substantially identical to that described in Ap- 
pendix A of [4"2"] . Results using this version of the code 
first appeared in [55] . The discretization method we use 
for the derivatives that appear in Eq.[6]is a standard one; 
e.g., for a field cj> at the grid location {i, j, k} on an x, y, z 
grid, we would have 

^xVx<t>i,j,k = h ~ 2 (<Pi+i,],k + 4>i-i,jM - %<l>i,j,k) (Al) 
1 . 



H+\,j+\,k — (Pi+l,j-l,k 
■<Pi-l,j+l t k + <f>i-l,j-l,k) (A2) 



Appendix B: Two point clustering statistics in real 
and redshift space 

Calculating the two point correlation function for TV 
particles by direct pair counting requires ./V 2 operations. 
Considering the large number of particles used in the 
simulation we make use of an estimator introduced by 
[7TI 172"] . In this approach a density field is constructed on 
iV grid cells and the correlation function is then calculated 
as 



£(|ryl) = 



N inn N p (\ nj \) ^ 



£$>( r <)*fa))* (Bl) 



and h = 2™, with n representing the level of refinement 
in the multigrid relaxation. 



where = (n(ri)— < n >)/ < n > is the density 

fluctuation in the i th bin of the grid. The sum extends 
over all N p pairs separated by distances between r— Ar/2 
and r+Ar/2. We also sum over eight realisations for each 
simulation of a particular cosmological model, iV run = 8. 
This procedure scales as N 2 lid and requires fewer opera- 
tions then direct pair counting as usually iV gr id <C N. 

This approach is a far more efficient method to mea- 
sure the correlation function than direct pair counting 
and has been shown to be extremely accurate at repro- 
ducing the full two point function for a range of grid sizes 
[see Figure 5 in [72] . [73] and [7J also used a grid based 
calculation with FFT to measure the correlation func- 
tion and found that this method accurately reproduces 
the £(r) found from direct pair counting. Using a lower 
resolution simulation, we have verified that the using the 
estimator in Eq. |B1| reproduces the correlation function 
measured using the standard |75j method. 

This grid based method limits the accuracy of our 
measurements to scales larger than a few grid cells, 
r > L box /N srid . For the L box = 400 Mpc//i simulation we 
use a iV grid = 200 3 grid while for the L box = 256 Mpc/h 
simulation we use N gIid = 256 3 . Errors on the z = 
correlation function represent the scatter amongst 8 re- 
alisations of the same cosmology where different random 
number seeds where used to generate the initial condi- 
tions for the simulations. The errors on the measure- 
ments at z > were obtained by Jackkife sampling from 
a single simulation by dividing the simulation volume 
into iV 3Ub = 8 equal subvolumes and then systematically 
omitting one subvolume at a time in order to calculate 
the correlation function on the remaining N snb — 1 volume 
[see [75] for more details of this method] . The redshift 
space correlation function is obtained from the simula- 
tions after averaging over the £(s) obtained by treating 
the x, y and z directions in turn as the lines of sight. 
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